rm(list=ls())
setwd("/Users/afanlo/Documents/Projects/Kroenig Rep")
library(dplyr)
library(MASS)
library(xtable)
library(stargazer)
library(lmtest)
library(sandwich)

## Upload Kroenig Data
cutoff_df<-read.csv("DF_WithCutoffs.csv")
cutoff_df<-cutoff_df[,c(3:9, 29:148)]
dataframe<-read.csv("data_may2020.csv")
df<-merge(dataframe, cutoff_df, by=c("ccode", "ccode2", "crisno", "year", "crisname",
                                     "abbrev1", "abbrev2"))

### Creating a function to compute clustered standard errors
vcovCluster <- function(
    model,
    cluster
)
{
  require(sandwich)
  require(lmtest)
  
  cluster = as.factor(cluster)
  
  if (!is.null(model$na.action)){
    omit.rows = model$na.action
    cluster = cluster[-omit.rows]
  }
  if (nrow(model.matrix(model)) != length(cluster)){
    stop("something's not working: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap or randomization inference instead).")
  }
  
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- na.omit(apply(estfun(model), 2, function(x) tapply(x, cluster, sum)))
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

###
### NUCLEAR RATIO WITH APPENDIX F - Fanlo/Sukin - RE-CODINGS 
##

## Changing back to ICB Coding of 1969 Sino-Soviet Border War
df$victory<-ifelse(df$crisno==231, 0, df$victory)

##Changing Berlin Wall Coding to a Stalemate
df$victory<-ifelse(df$crisno==185,  0, df$victory)

## Don't need to change Korean War outcome, since the starting point is Kroenig's data and he
## already recoded that outcome as a stalemate

## Ratio Table
table(df$nucratio)


## Models Begin
mod1<-glm(victory~nuc_rat + proximity + polity21 + 
            capshare + strikea +  tpop_1 + viol + cris_ave, data=df, family=binomial(link="logit"))
mod1test<-coeftest(mod1, vcov=vcovCluster(model=mod1, cluster=factor(df$crisisdyad)))

###
## NUCLEAR RATIO WITH ORIGINAL ICB OUTCOME Codings
###

## Upload Kroenig Data
cutoff_df<-read.csv("DF_WithCutoffs.csv")
cutoff_df<-cutoff_df[,c(3:9, 29:148)]
dataframe<-read.csv("data_may2020.csv")
df<-merge(dataframe, cutoff_df, by=c("ccode", "ccode2", "crisno", "year", "crisname",
                                     "abbrev1", "abbrev2"))

## Changing back to ICB Coding of 1969 Sino-Soviet Border War
df$victory<-ifelse(df$crisno==231, 0, df$victory)

##Changing back to ICB coding of Korean War
df$victory<-ifelse(df$crisno==133 & df$ccode==365,  1, df$victory)

## Models Begin
mod2<-glm(victory~nuc_rat + proximity + polity21 + 
            capshare + strikea +  tpop_1 + viol + cris_ave, data=df, family=binomial(link="logit"))
mod2test<-coeftest(mod2, vcov=vcovCluster(model=mod2, cluster=factor(df$crisisdyad)))

###
## NUCLEAR RATIO WITH KROENIG'S OUTCOME CODINGS
###

## Upload Kroenig Data
cutoff_df<-read.csv("DF_WithCutoffs.csv")
cutoff_df<-cutoff_df[,c(3:9, 29:148)]
dataframe<-read.csv("data_may2020.csv")
df<-merge(dataframe, cutoff_df, by=c("ccode", "ccode2", "crisno", "year", "crisname",
                                     "abbrev1", "abbrev2"))

## Models Begin
mod3<-glm(victory~nuc_rat + proximity + polity21 + 
            capshare + strikea +  tpop_1 + viol + cris_ave, data=df, family=binomial(link="logit"))
mod3test<-coeftest(mod3, vcov=vcovCluster(model=mod3, cluster=factor(df$crisisdyad)))



## Create Table 1 from Appendix H 
library(stargazer)
stargazer(mod1test, mod2test, mod3test)

